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RESEARCH  SUMMARY 


Since  October  I,  1973,  we  have  been  supported  by  ARPA  IPT  office  to  carry 
out  research  in  image  analysis  and  modeling.  The  emphasis  of  our  research  in 
the  past  two  years  has  been  on  the  development  of  suitable  mathematical  models 
for  images  which  are  useful  for  image  processing  tasks  such  as  efficient  cod- 
ing, enhancement,  recognition,  and  information  extraction.  Our  accomplish- 
ments have  been  recorded  in  detail  in  our  progress  reports.  We  summarize 
here  some  of  the  highlights. 

A.  Preprocess i ng 

We  have  obtained  significant  results  in  both  image  restoration  and  ef- 
ficient coding. 

1 . Image  phase  [l] . 

Although  it  is  well  known  that  the  phase  of  the  Fourier  transform  of  an 
image  is  generally  more  important  than  the  magnitude,,  most  past  work  on  two- 
dimensional  digital  filter  design  concentrated  on  magnitude  specification  only. 
We  demonstrated  that  the  phase  accuracy  of  image  processing  filters  are  ex- 
tremely important.  Even  if  the  desired  filter  has  linear  phase,  failure  to 
specify  it  may  lead  to  disaster.  We  also  developed  methods  of  designing  two- 
dimensional  digital  filters  which  specify  both  the  magnitude  and  the  phase. 

2.  Recursive  estimation. 

By  modeling  images  as  two-dimensional  random  fields,  it  becomes  feasible 
to  apply  the  Kalman  formulation  of  recursive  estimation  to  image  restoration. 
However,  in  trying  to  derive  the  optimum  two-dimensional  estimator,  one  en- 
counters fundamental  mathematical  difficulties.  We  have  nevertheless  develop- 
ed several  suboptimum  estimators  which  in  practice  perform  almost  as  well  as 
the  optimum. 
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3.  Iterative  image  restoration  [2,3,4]. 

Many  image  degradations  can  be  approximated  by  linear  model;,.  Then 
image  restoration  on  the  computer  becomes  the  problem  of  solving  a set  of 
linear  algebraic  equations.  Because  of  Image  noise,  conventional  solution 
methods  are  unusable.  We  have  developed  an  iterative  method  (the  projection 
method)  for  doing  the  restoration  which  offers  a tradeoff  between  noise  and 
image  sharpness.  Compared  to  the  singular  value  decomposition  method  studied 
at  Purdue  and  elsewhere,  this  Iterative  method  offers  a similar  performance 
with  drastically  reduced  computation  requirements. 

4.  Three-dimensional  reconstruction  [53. 

We  have  studied  several  methods  of  reconstructing  three-dimensional 
structures  from  two-dimensional  x-ray  pictures.  Specifically,  we  have  in- 
vestigated the  effects  of  quantization,  beam  divergence,  and  unknown  beam 
strength.  It  was  found  that  64  levels  of  quantization  Is  sufficient,  that  a 
beam  divergence  of  up  to  10°  can  be  tolerated,  and  that  an  unknown  beam 
strength  introduces  a ring-1  Ike  structure  to  the  reconstruction. 

5.  Error-free  DPCM  codes  for  ERTS  Imagery  [6,7]. 

The  difference  statistics  of  ERTS  imagery  have  been  measured;  and  based 
on  these  several  easily  implementable  classes  of  DPCM  codes  have  been  develop- 
ed. They  reduce  the  bit  rate  from  8 bits  per  picture  element  to  about  4 bits 
per  picture  element.  A slightly  more  complicated  adaptive  code  reduces  the 
bit  rate  to  about  2.5  bits  per  picture  element.  These  codes  are  error-free  in 
the  sense  that  they  do  not  Introduce  any  distortion  to  the  images. 

B.  Image  Segmentation 

We  have  taken  two  approaches  to  image  segmentation:  edge  extraction,  and 

region  growing. 
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I.  Edge  extract  ion  [8], 


We  have  developed  several  Improved  local  operators  for  edpe  extraction. 
They  work  quite  well  even  In  a noisy  and  blurred  environment.  These  operators 

have  been  applied  to  change  detection  of  missile  imagery  with  excellent  results, 

2.  Region  growing  [9]. 

We  have  developed  a computer  algorithm  called  BLOB  which  segments  an  image 
•nto  regions  so  that  points  in  the  same  region  have  similar  characteristics. 

The  BLOB  algorithm  was  used  to  Increase  the  accuracy  of  classifying  multi- 
spectral  ERTS  data.  These  data  had  been  classified  point  by  point  using 
spectral  signature.  By  the  application  of  BLOB,  regional  classification  be- 
came possible.  This  not  only  increased  the  classification  accuracy  (by  about 

5?  on  the  average)  but  also  reduced  the  classification  time  (by  a factor  of 
of  about  30:1  ). 

3.  Image  decomposition  [8], 

Inage  decomposition  can  be  considered  as  a generalization  of  image  seg- 
mentation. In  the  latter,  an  image  is  segmented  into  non-overlapping  regions: 
while  in  the  former,  an  image  is  decomposed  into  a sum  of  generally  over- 
lapping images.  For  example,  an  image  might  be  decomposed  into  the  sum  of  a 
low-spatial-frequency  image  and  a h i gh-spa ti a I- frequency  image.  Many  image 
processing  tasks  are  facilitated  if  the  image  is  decomposed  into  simpler  com- 
ponents and  each  component  is  handled  according  to  its  own  characteristics. 

We  have  been  developing  algorithms  to  decompose  an  image  into  three  components: 
edges,  background  ( slowly-changing),  and  textures.  This  idea  of  image  de- 
composition has  been  applied  to  image  noise  reduction.  Applying  a Wiener 
filter  to  a noisy  image  reduces  the  noise  but  also  blurs  the  edges  of  the 
objects  in  the  image.  By  treating  the  edges  separately,  one  is  able  to  re- 
duce the  noise  and  in  the  same  time  retain  edge  sharpness. 
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C.  Pattern  Classification 


1.  Texture  analysis  [10]. 

We  have  developed  a set  of  image  texture  descriptors  based  on  extrema 
along  scan  lines.  These  texture  descriptors  are  very  easy  to  measure,  yet 
perform  as  well  as  or  even  better  than  conventional  (tedious  to  calculate) 
descriptors  (such  as  the  spatial-dependence  matrix)  in  most  pattern  recogni- 
tion tasks. 

2.  Optimum  feature  subset  selection  algorithms  [11]. 

Patterns  to  be  classified  by  a recognition  system  are  generally  character- 
ized by  a set  of  measurements  or  features.  Often,  the  dimensionality  of  this 
feature  space  is  too  large  for  efficient  and  reliable  application  of  existing 
classification  techniques.  The  feature  extraction  problem  is  then,  to  reduce 
the  dimensionality  of  the  feature  space,  without  signif icfintly  affecting  the 
discriminatory  capacity  cf  the  feature  set.  One  approach  to  feature  extraction 
is  tj  select  a smaller  subset  of  m features  of  the  set  of  n original  features 
(m<n).  Exhaustive  search  is  computationally  unfeasible  even  for  modest  values 
of  n and  m.  Using  a branch  and  bound  approach,  we  have  developed  algorithms 
which  are  efficient  and  guaranteed  to  be  optimal.  These  algorithms  were 
applied  to  choose  the  best  subset  of  12  out  of  24  channels.  There  are 
2,704,150  candidate  solutions.  The  algorithms  obtained  the  best  subset  with 
the  computational  effort  equivalent  to  computing  the  criterion  for  6000  sub- 
sets. The  savings  are  indeed  substantial. 
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IMAGE  DECOMPOSITION 


T.  S.  Huang  and  J.  W.  Burnett 

I . Introduction 

In  previous  reports  [1],  [2]  we  have  modeled  image  scan  lines  as  Markov 
jump  processes.  This  model  led  to  nonlinear  noise  reduction  and  image  seg- 
mentation algorithms  that  are  superior  to  linear  techniques  currently  in  use. 

The  recursive  calculation  of  a conditional  probability  involving  the 
boundary  component  of  the  scan  line  was  the  key  to  the  nonlinear  algorithms. 
Once  this  conditional  probability  had  been  calculated  the  scan  line  could  be 
segmented  using  a maximum  likelihood  approach  or  noise  could  be  reduced  with 
either  maximum  liklihood  or  minimum  mean  square  error  estimators.  Further, 
since  the  conditional  probability  obeyed  a recursive  relationship  the 
execution  time  and  memory  requirements  of  these  algorithms  were  kept  to  a 
mini  mum. 

In  this  report  we  start  to  extend  the  one  dimensional  results  to  two 
dimensions.  The  first  step  in  this  extension  is  the  recursive  estimation  of 
a two  dimensional  constant. 

I I . Derivation 

Assume  we  have  the  observations 
(1)  R(x,y)  ■=  Ixy  + BW(x,y) 

where  I is  one  of  K possible  values  a]f  a2,  ...,  aK  and  W(x,y)  is  a two 
dimensional  Wiener  process  (the  integral  of  white  Gaussian  noise.) 

Equation  (1)  is  sometimes  written  as 

O')  dR  = Idxdy  + BdW  or 
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(,,,)  1 = 1 + BN(x.y) 

where  N(x,y)  is  white  Gaussian  noise. 

Let  S = 12.  t =31 
rn  n qn  n 


r,q  = 1,  2,  3 n 


^ 1 » 2,  3,  ... 


Def i ne 


^qn  " R(Sm’tqn)  + R ( 5 r- 1 n * tq - 1 n } ' R(Srn*tq-ln)  “ ^r-ln’^J 


From  (1) 


n = + BAW 

rqn  ^2  rqn 


r.q  = 1,  2 


where 


Awrqn  = W(Srn’tqn)  " ^r-ln'V  ~ W($rn  ’ V 1 n}  + W(Sr-|n*tq-ln) 

Now  AW^  is  a linear  combination  of  Gaussian  random  variables  and  hence 
Gaussian.  Further,  since  W is  an  independent  area  process  AW^  is  indepen- 
dent of  AWgbn  for  a^r , b/q.  Thus 


'rqn  2 
n 


^*9“  1*  2,  • • * 9 n 


are  independent  zero  mean  Gaussian  random  variables  with 

Let  = P ( I = a . | R(S  ,t  ),  r,q  = 1,2  n 

J J 1 rn’  qn  * n 


var i ance 


= conditional  probability  the  random  field  I has  value  a 

j 

given  the  observations  R upto  the  point  x,y. 

Note  that  nrqn  is  obtained  from  R by  an  operation  that  (except  for  some 
boundary  conditions)  is  nonsingular.  Thus  if  we  define  R(x,y)  = 0 for  x < 0 
or  y < 0 then 


\ 


p ( I " aJ  |R(S  »t  ))  = P(I  n a.  In  ) 
rn  qn  j ' re" 


From  Bayes  rule  we  have 


(n) 

P:(x,y)  = 


P.  (0,0)  exp [ j 

J a r\£- 


n 

z (n 


2B  xy  r,q=l 


rqn 


ojxy  2 


2 P£(0,0)  exp [ Z (n  - !£I)2] 

M 2 J2xy  r#q-l  rqn  IF 


(2) 


where  P ^ (0 , 0)  is  the  a-priori  probability  that  I = aj. 
Equation  (2)  can  be;  rewritten  as 


P.(0,0)  exp [-4-  Z n a. 
(n)  J B2  r,q=l  rqn  J 

P (x,y)  

J i •*» 

z.  P«  (0,0)  exp[~  Z n a0 
i B2  r,q-l  rqn  * 


2 

vy] 


2B 


4*V] 


2B 


By  definition  of  the  two  dimensional  stochastic  integral  [3], 
1 a.  n a’ 


im 


J s n 


n-*00  B r,q=l 


rqn 


fo  fo  dR(u»v) 


-i-  R(x,y) 
B 


Sine*1  exp  (•)  is  a continuous  function 


2 

a xy 


P. (0,0)  exp[B  2a.R(x,y)  ~bri 
1 im  P.  (x,y)  - P . (x,y)  ?B 

n-KJO  J J ^ 


n-K» 


aT  ajxy 

Z Pj  (0,0)  exp  [— »■  R(x,y) £*— 

JM  B*  2B2 


] 


Each  P is  of  the  form 


fj  (x,y) 

^TxTyT 


where  f£  ■»  P£(0,0)  exp  [B~2a£R(x,y)  - 


a£xy 


2B 
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Using  a result  of  Wong  [3],  [4]  it  can  be  shown  that  each  f.  obeys  a recur 
sive  relationship 


df.(x,y)  = B 1 a . f . (x,y)dR(x,y)  + 


B'V’Mx.y )/V  dRCf^yJdRfx,^) 


(b) 


Thus  an  incremental  change  in  each  f ^ can  be  calculated  by  knowing  the 
"present  value"  of  f.  and  the  present  value  of  the  observation  (see  Fig.  1). 

III.  Conclusions 


Equation  (4)  is  significant  in  that  it  gives  a two  dimensional  recursive 
method  for  calculating  a conditional  probability.  From  (4)  we  can  recur- 
sively calculate  each  f^  and  hence  each  . While  this  result  may  not  find 
widespread  application  it  nonetheless  represents  the  first  step  in  deriving 
a recursive  expression  for  a conditional  probability  when  the  random  field  I 

(1)  is  a Markov  jump  process. 
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Figure  I Recursive  Calculation  of  df(x,y)  from 

past  values  of  f(x,y)  and  present  values 
of  dR(x,y) 
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DIGITAL  EDGE  RESTORATION  IN 
LINEARLY  FILTERED  IMAGES 

D.  P.  Panda  and  A.  C.  Kak 

I . INTRODUCTION 

Linear  minimum  means  square  error  (MMSE)  estimation  techniques  have  com- 
monly been  used  in  restoring  images  degraded  by  the  imaging  system  and  observed 
in  the  presence  of  additive  noise  (see  [1]  for  brief  survey).  These  techniques 
often  distort  the  "edges'1  of  the  restored  image.  The  cause  of  the  distortion 
is  the  fact  that  the  edges  mark  the  boundaries  between  two  adjacent  "clusters" 
of  different  local  statistical  properties  and  the  estimation  of  a picture 
element  (pixel)  in  one  cluster  nearby.  Depending  on  the  intended  use  of  the 
restored  image  the  edge  distortion  may  or  may  not  be  too  undesirable.  If  the 
image  is  meant  for  the  human  visual  system  (HVS)  then  it  is  desirable  to  have 
sharp  and  undistorted  edges  in  the  image.  This  is  so  because,  as  Hunt  [2] 
pointed  out,  the  frequency  response  of  the  HVS  is  such  that  it  differentiates 
the  lower  spatial  frequencies  and  amplifies  the  higher  spatial  frequencies 
which  amounts  to  emphasizing  the  edge  contents  of  the  image.  When  a linear 
MMSE  restoration  filter  does  not  have  any  edge  preserving  constraint  or  measure 
it  shall  be  necessary  to  find  a post-processing  technique  which  would  enhance 
the  edges  in  the  restored  image  (the  term  restored  image  is  used  here  to  mean 
the  output  of  the  linear  MMSE  restoration  filter).  Many  frequency  domain 
methods  are  already  available  for  edge  enhancement.  One  such  frequency  domain 
technique  is  due  to  Schreiber  [3]  in  which  he  subtracted  from  the  image  a 
weighted  low-pass  version  of  it  to  increase  the  relative  amplitudes  of  the  high 
spatial  frequency  components  and  restore  some  of  the  sharpness.  Frequency 
domain  operations  may  be  preferred  when  being  implemented  optically.  But  in 
digital  implementation  spatial  domain  operations  are  expected  to  be  more 
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attractive  than  frequency  domain  operation  from  the  storage  and  computation 
point  of  view.  We  present  here  a spatial  domain  post-processcr  for  the  en- 
hancement of  the  edges  in  the  images  restored  linearly  by  the  MMSE  criterion. 

1 1 * EDGE  DISTORTION  IN  LINEARLY  FILTERED  IMAGES 

There  are  several  ways  the  edges  in  an  image  can  be  modelled.  Each  of  the 
many  edge  detection  schemes  devised  by  various  researchers  is  based  implicitly 
on  a specific  model  of  the  edges.  In  terms  of  power  spectral  density  the 
edges  contribute  to  the  relatively  higher,  spatial  frequency  components  of  an 
image.  In  some  spatial  domain  edge  detection  methods  an  edge  is  defined  as  the 
set  of  pixels  that  mark  the  transition  of  the  local  mean  grey  level  from  one 
value  to  another,  while  in  other  methods  the  change  in  the  variance  of  the  grey 
levels  is  also  taken  into  account.  Whether  it  is  the  first,  the  second,  or  the 
higher  order  statistics  that  are  used  for  edge  detection  it  is  obvious  that  the 
edges  are  the  boundaries  of  regions  that  contain  different  local  statistical 
properties.  When  these  local  properties  are  used  as  features  in  a clustering 
procedure  the  above  mentioned  regions  cluster  together  in  the  feature  space. 
These  clusters  in  the  feature  space  are  separated  from  one  another  by  the  edges 
in  the  image  space  and,  in  general,  correspond  to  different  objects  in  the  ob- 
ject space  (see  Fig.  1).  Researchers  in  the  area  of  remote-sensing  make  use  of 
this  concept  in  the  spatial  clustering  of  multi-image  data  and  scene  classifi- 
cation (see,  for  example,  Haralick  and  Dinstein  [5]  and  Gupta  and  Wintz  [^]). 

We  shall  use  here  the  terms  edges,  "cluster  boundaries",  and  boundaries  inter- 
changeably. 

When  a linear  estimator  is  used  for  filtering,  the  local  properties  of  an 
estimated  pixel  are  affected  by  the  local  properties  of  all  the  pixels  included 
in  the  observation  set  for  estimating  that  pixel.  Therefore,  in  the  neighbor- 
hood of  the  boundaries  of  two  clusters  the  estimated  pixels  will  have  local 


13 


properties  depending  on  the  local  properties  of  both  the  clusters  rather  than 
just  the  cluster  to  which  the  pixel  belongs. 

Def ini tions 


Let  s(p,q)  be  the  grey  level  of  the  Image  at  the  point  (p,q).  Suppose 
the  image  is  scanned  in  certain  direction  6,  0 < 6 < n.  In  any  particular 
scan  line  we  define  a run  R as  the  set  of  consecutive  pixels  Such  that  the 

grey  level  of  those  pixels  is  monotomical ly  increasing  or  decreasing.  If 

(PpV  and  (p2’q2)  are  the  two  end  P°,nts  of  the  run  then  we  define  the  con- 
trast of  the  run  as 


C(R)  = Mpj.qj)  - s(p2,q2)|/[s(p),q1)  + s (p2,q2)J 
and  for  any  point  (p,q)  in  the  run  we  define  sharpness 


^(p.q) 


C(R) 

((P,-P2)2  + (q,-q2)2),/2 


We  define  the  sharpness  error  criterion  between  an  image  and  its  estimate  as 
e(p,q)  = l (p,q)  - Y (p,q)]2 

e=e|  0 0 

where  * indicates  the  value  in  the  estimated  image,  6| , 02 are  k dlf 

ferent  directions  along  which  the  images  are  to  be  scanned  for  computing  the 
error. 


Ml.  EDGE  RESTORATION 

Since  the  effect  of  pixels  from  one  cluster  on  the  pixels  of  a neighboring 
cluster  near  the  boundary  Is  the  cause  of  the  edge  distortion,  If  we  introduce 
in  the  linear  filter  nonlinearity  neat'  the  boundary  such  that  an  estimated 
Pixel  of  a cluster  Is  a linear  combination  of  pixels  of  only  that  cluster  then 
we  should  expect  a great  reduction  in  the  distortion.  But  if  no  such  attempt  Is 
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made  to  prevent  the  edge  distortion  In  the  linear  filter  then  a post-processor 
is  needed  to  correct  the  possible  distortion  In  the  image.  The  post-processor 
presented  here  does  this  correction  in  the  following  steps.  First  it  detects 
in  the  output  of  the  linear  filter  the  boundaries  between  adjacent  clusters. 
Then  it  changes  the  grey  levels  and  the  gradients  of  the  pixels  in  those  re- 
gions in  such  a fashion  that  the  sharpness  of  the  boundary  is  increased.  If 

needed,  it  smooths  out  further  the  regions  that  do  not  constitute  the  boundary. 
Edge  Detection 

Let  s (p j , q j ) and  s(p2>q2)  be  the  grey  levels  at  pixels  (p] , q ] ) and  (p  q ) 

in  adjacent  clusters  1 and  2,  respectively.  We  will  assume  that  s (p  ^ , q ^ ) and 

s(p2,q2),  respectively,  are  Gaussian  with  means  y ^ and  y2  and  variances  o2  and 
2 

a2'  A,so»  f°r  tractabl 1 i ty  In  deriving  the  boundary  detection  algorithm  we 
will  replace  the  expression  for  contrast  of  a run  R by  the  following 

local  contrast  C(R)  = (s(k,A)  - s(i,j))/2y  (1) 

where  y is  the  average  grey  level  of  the  region,  and  (k,Jt)  and  (i,j)  are  the 
maximum  and  the  minimum  end  points,  respectively,  of  the  run  R.  When  the  run 
R is  across  a boundary  those  points  will  not  be  in  the  same  cluster.  Let 

D(R)  = s (k, £)  - s (i , j ) 

For  a run  R across  a boundary  the  conditional  probability  density  of  D given 
that  a boundary  has  been  detected  is  given  by 

PD(d|B)  " N(Mb,ab2)  (2) 

where  B is  the  event  that  an  edge  has  been  detected,  y and  a 2 are  the  mean 

b b 

and  the  variance  respectively  of  the  Gaussian  density  function.  These  para- 
meters are  given  by 
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H,  = m2  - yi 


2 2 2 
Gb  = ai  + a2 

If  instead  the  pixels  (k,£)  and  (i,j)  belonged  to  the  same  cluster,  say  cluster 
1,  i.e.,  (k,£)  and  (i,j)  are  two  pixels  at  the  opposite  ends  of  a run  of  nono- 
tonically  increasing  or  decreasing  grey  levels  but  the  region  between  them  is 
not  a cluster  boundary,  then 

PD(d|C)  « N(uc,ac2)  {k) 

where  C denotes  the  event  that  the  run  is  inside  a cluster.  The  parameters  of 

2 

the  Gaussian  density  function  N(Mc,0c  ) are  given  by 


uc  = o 


2 . 2 
o = 2o, 

c 1 


If  it  is  assumed  that  in  the  output  of  the  linear  filter  the  runs  are  of 
approximately  equal  length  then  given  a particular  value  of  D(R)  whether  or 
rot  the  run  R is  across  a cluster  boundary  is  now  a pattern  recognition  prob- 
lem. If  P(B)  and  P (C ) , respectively,  are  the  prior  probabilities  of  the 
existence  of  a boundary  and  a cluster  In  a particular  region  then  the  Bayesian 
solution  to  the  pattern  recognition  problem  is 


f»>' 


imp! ies  B 


* 1 impl ies  B or  C 


< 1 implies  C 


If 


Using  appropriate  Gaussian  densities  for  p(d|B)  and  p(d|c)  we  get  the  dis- 
criminant functions, 


(d)  = ' [V  JT  “ J ("T  + *og  o?)  + log  P (B) 

% % % 


, d^  dp  . p 

gc(d)  "t  j + ^r"  ^t  + 109  °c)  + ,og  p(c)  (7) 

C C C 

In  terms  of  the  discriminant  function,  (6)  can  be  expressed  as 
9^ (d)  - gc (d)  > 0 imp] i es  edge 

< 0 implies  no  edge  (8) 

There  exists  a tradeoff  between  the  accuracy  of  the  detection  of  boundary 
and  the  speed  and  the  cost  of  computation  involved.  Where  as  using  both  the 
first  and  the  second  order  statistics  is  expected  to  give  better  performance 
cost.  Often  just  the  first  order  statistics  gives  satisfactory  results  with 
less  cost.  Hence,  we  take  the  variance  02  and  o2  to  be  equal  (=  o2)  thus  re- 
ducing the  complexity  of  the  discriminant  function  as  shown  below 

2 2 . 2 

0.  = o = 2o 
b c 


Therefore 


%<d> 


Vb) 


dMb  Mh 

— 2 J + log  p(b) 

2 02  402 


duc  Uc 

— J ' —T  + J°9  p(c) 
2o  in 


The  decision  rule  now  is 


d > t *►  edge 


< t -*■  no  edge 


where  the  threshold  t is 


t 


2 2 

V^c 

4o2 


+ log 


P(c) 

HbT 


2o 

The  probability  of  error  of 
little  computational  complexity, 


the  classifier  is  decreased,  at  the  cost  of 
by  changing  the  decision  rule  as  follows; 


a 


9g(d)  - gc  (d)  > t^->  boundary 


< t no  boundary 
otherwise  undecided 


When  undecided  a new  run  of  monotonical ly  increasing  or  decreasing  grey  levels 
is  added  to  the  previous  run.  This  process  is  continued  till  a boundary  is 
detected.  The  width  of  the  boundary  is  the  cumulative  length  of  all  runs  used 
in  making  the  decision.  Since  the  edge  width  in  the  output  of  a linear  filter 
does  not  distort  without  limit  an  empirically  determined  threshold  can  be  set 
for  the  cumulative  run  length.  If  at  any  time  in  a decision  making  the  cum- 
ulative run  length  meets  the  threshold  "no  boundary"  is  decided. 

Boundary  Enhancement  and  Cluster  Smoothing 

The  sharpness  of  the  edge  can  be  Increased  either  by  increasing  local  con- 
trast or  by  reducing  the  boundary  width.  To  increase  the  contrast  by  a small 
fraction  we  must  increase  or  decrease  the  grey  levels  in  the  boundary  region 
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to  a great  extent  which  will  be  "out  of  step"  with  the  rest  of  the  pixels  in 
the  ciuster,  thereby  creating  artifacts.  For  example,  if 


s2(M)  = 75 
= 25 

contrast  = = i/2 


Increasing  s (k,£)  5 times  to  375  r.. 


decreasing  s f ( i , j ) 5 times  to  5 gives 


or 


contrast 


= 375-25  = 350 
375+25  kOO 


7/8 


contrast  7/8 

which  is  even  less  than  a twofold  improvement.  But  decreasing  the  boundary 
width  to,  say,  1 gives  a very  large  order  of  improvement  in  the  sharpness. 

This  improvement  is  proportional  to  run  length.  The  worse  the  edge  distortion 
is  the  larger  is  the  run  length  and  consequently  the  higher  is  the  improvement 
in  the  sharpness  of  the  edge. 

The  boundary  width  is  reduced  by  a nonlinear  mapping  of  the  grey  levels  by 
which  the  grey  levels  above  the  midpoint  get  mapped  to  s2(k,£)  and  those  below 
the  midpoint  get  mapped  to  Sj(i,j). 

The  next  step  after  enhancement  of  the  edges  is  the  smoothing  of  the  in- 
terior of  the  clusters.  This  smoothing  may  or  may  not  be  necessary  depending 
on  the  residual  noise  present  in  the  iineariy  filtered  image.  If,  for  example, 
the  linear  filter,  such  as  the  one  by  Habibi  [6],  is  the  kind  that  used  the 
image  correlation  coefficients  as  its  parameters  then  the  amount  of  residual 
noise  in  a ciuster  in  the  filtered  image  might  depend  on  how  well  the  corre- 
lation coefficients  used  in  the  filter  design  match  that  of  the  cluster  under 


19 


consideration.  Usually  the  come, at, on  coefficients  have  to  he  estimated  fro. 
one  or  a few  sample  images  and  unless  the  restoration  filter  Is  adaptive  It 

b°  ,mPOSSib'e  f°r  ,l  “ ‘I-  c-usters  encountered.  «hen  needed, 

the  smoothing  Is  done  In  the  post-processor  by  a simple  first  order  digital 

recursive  filter.  The  pole-zero  locations  of  the  filter  may  be  controlled  to 
get  the  desired  amount  of  smoothing. 
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Figure  k The  image  of  Fig.  3 
restored  along  the  horizontal 


after  the  edges  have  been 
and  the  vertical  directions 


Figure  5 The  image  of  Fig.  l>,  with 
probability  mapped  to  nearest  grey 
probabi 1 i ty 
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IMAGE  SEGMENTATION  BY  UNSUPERVISED  CLUSTERING  III 
M.  Y.  Yoo  and  T.  S.  Huang 

1.  Introduction 

This  is  a continuation  cf  the  work  reported  ir,  the  last  two 

interim  reports.  The  basic  philosophy  of  this  work  was  presented  in  detail 

before  and  only  the  new  clustering  algorithm  which  was  used  for  this  work  will 
be  discussed. 

2.  Clustering  Algorithm 

We  use  the  non-parametri c clustering  algorithm  developed  by  Fukunaga  et. 

al.  [1].  For  large  sample  sizes  this  algorithm  becomes  basically  the  same  as 
the  valley  seeking  algorithm  [2], 

The  only  difference  in  our  implementation  is  that  the  conditional  proba- 
bility density  function  (this  is  not  an  actually  normalized  probability 
function » the  value  of  the  function  denotes  the  number  of  points  which  have 
the  same  features  in  the  original  image)  in  the  feature  plane  is  available 
rather  than  sampled  vactors.  The  algorithm  consists  of  three  major  steps: 
i)  determination  of  neighborhood;  il)  searching  of  roots;  iii)  assignment  of 
parent. 

Let  X be  the  integer  pair  in  the  feature  plane  (our  features  are  quantiz- 
ed to  integer  levels),  and  D ()T)  denote  the  number  of  points  In  the  original 
picture  which  have  the  feature  X.  The  neighborhood  of  X is  the  set  defined  as 

ne(x)  - { Y | d(X,Y>  10,X  t 7} 

where  0 is  a preset  parameter  and  d means  Euclidean  distance. 

Then  the  total  number  of  the  original  image  points  which  has  the  feature 

in  the  neighborhood  n (x)  is 

0 


D(Y) 


Nfl  (X)  - Z 

YGnQ(jT) 

We  define  the  factor  G.  (X,Y)  as 

0 


% (X,Y)  = 


Ne (7)  - Ne  ^ 

dor.Y) 


Actually  Ge(X,Y)  Is  the  discrete  version  of  the  directional  derivative.  Then 
the  algorithm  is  follows.  If  n0  00  is  empty,  we  call  X root,  and  if  the 

neighborhood  of  X is  not  empty,  we  calculate  S (X",Y)  defined  as 

0 

S0 (X,Y)  = MAX  G0(X,Y) 

Yen0  ()0 


,f  Se(X,Y)  is  negative,  X is  also  called  root  and  if  Se(X,Y)  is  positive  we 
call  Y the  parent  of  X.  Sq(X,y)  = 0 case  Is  a little  bit  complex  and  we 
recommend  the  original  paper  for  details. 

Therefore,  every  point  in  the  feature  plane  are  either  classified  as 
roots  or  given  the  parent  relationship.  (roots  do  not  have  parents.)  Every 
point  except  roots  are  merged  into  one  of  the  searched  roots  by  parent  re- 
lationship and  the  number  of  roots  determined  that  of  clusters. 

3.  Experimental  Results 

The  number  of  clusters  is  very  much  dependent  upon  the  size  of  the  para- 
meter 6 which  is  used  when  we  define  the  neighborhood.  Smaller  6 gives  too 
many  artificial  clusters  and  larger  fl  classifies  the  whole  data  as  a cluster. 

The  appropriate  size  of  6 for  this  work  was  3.  For  large  size  data 

there  is  a serious  storage  problem,  and  some  preprocess  for  "throwing  away"  of 

irrelevant  data  is  highly  necessary.  We  shall  present  some  experimental  results 
in  a late'*  report. 
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TEXTURE  BOUNDARY  DETECTION 
0.  R.  Mitchell  and  W.  K.  Chan 

In  an  attempt  to  develop  efficient  image  segmentation  algorithms,  we  are 
concentrating  on  texture  boundary  detection.  Most  edge  detection  algorithms  do 
not  use  texture  information  at  all.  Our  initial  algorithms  use  texture  infor- 
mation exclusively.  The  combination  of  texture  and  gradient-type  boundary 
detectors  will  be  done  following  our  development  of  the  texture  edge  detection. 

The  present  procedure  uses  the  Max-Min  method  of  texture  analys is  descr i bed 
elsewhere  in  this  report.  However,  instead  of  using  the  method  for  a texture 
classification,  a max-min  feature  computation  is  performed  in  two  opposite  di- 
rections from  a point  in  the  picture.  The  resulting  features  are  compared.  If 

the  features  in  opposite  directions  differ  significantly,  we  mark  the  center 
point  as  a candidate  for  a texture  edge.  This  calculation  is  performed  itera- 
tively as  we  move  the  center  point  across  the  picture.  Peaks  in  feature 
differences  indicate  texture  edges  as  shown  In  Fig.  1.  More  resolution  to  the 

edges  is  found  by  repeating  the  scan  in  directions  perpendicular  to  suspected 
edges. 

Future  improvement  to  be  added  are  the  automated  selection  of  scan  length, 
L,  based  on  a regional  autocorrelation  and  the  syntactic  joining  of  loosely 
connected  edge  points  to  define  closed  boundaries. 
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MAX-MIN  MEASURE  FOR  IMAGE  TEXTURE  ANALYSIS 
0.  R.  Mitchell  and  W.  A.  Boyne 


I.  BASIC  MAX-MIN  MEASURE 

"his  project  is  a continuation  of  that  described  In  the  May  l-July  31, 
1975  interin.  report  [I],  This  texture  measure  is  based  on  the  intuition  that 
important  texture  information  Is  contained  in  the  relative  frequency  of  local 
extremes  in  intensity.  The  grey  levels  values  along  a one-dimensional  scan 
direction  are  first  sent  through  a smoothing  process  which  eliminates  reversal 
of  small  amplitude,  thereby  retaining  only  the  principal  extrema. 

The  smoothing  algorithm  is  described  as  follows;  let  xfc  be  the  grey  level 
of  the  kth  point  along  the  scan  line  and  let  V|<  be  the  ’'smoothed"  value.  Le, 

T be  the  value  of  a preassigned  threshold  parameter.  Let  us  start  with  y -x 
and  proceed  according  to  the  algorithm  shown  below: 


IF 

THEN 

yk  < Vi  - 7 

Vi  = Vl  ~ 2 

T „ T 

Vi  T<yk<  Vi  *2 

V>  = yk 

Vi  + 7 < vk 

yk+f  = Vl  + J 

figure  i illustrates  the  way  In  which  the  smoothing  process  eliminates  rever- 
sals of  .nail  amplitude.  Implementation  of  the  extrema  detecting  algorithm  Is 
shown  in  ,he  flow  chart  in  Fig.  2.  B?  repeating  this  process  for  several 
threshold  settings,  a group  of  extrema  counts  can  be  obtained  to  characterize 
the  texture.  An  example  Is  shown  In  Fig,  3 using  three  different  threshold 
levels.  The  set  of  numbers  (6, 10, Ml)  would  characterize  this  "texture". 

To  make  the  features  Invariant  to  multiplicative  (gain)  changes  the  logs- 
r i thm  of  the  data  is  used.  Since 
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log  a - log  b = log  k - log  k, 

a b 

the  extrema  count  in  Fig.  3 is  unchanged  following  a gain  change. 

I I . CLASSIFICATION  RESULTS 

The  texture  data  used  has  been  described  previously  [1],  Forty-nine 
64x64  samples  of  each  of  eight  textures  were  used.  The  ratios  of  the  number  of 
extrema  at  each  selected  threshold  to  that  at  other  thresholds  were  used  as 
features.  The  threshold  settings  were  chosen  empirically  to  be  130,  110,  90, 
70,  50,  30,  and  10.  (The  data  was  already  in  log  form,  ranging  from  0 to  255.) 
Training  set  data  showing  the  feature  values  for  a few  of  the  texture  samples 
i s shown  in  Tabl e 1 . 

The  49  samples  of  each  texture  were  divided  into  36  training  samples  and 
13  test  samples.  Then  six  features  were  calculated  for  each  sample  and  two 
classification  techniques  were  used:  (l)  a simple  normalized  Euclidean  dis- 

tance measure  with  all  features  weighted  equally,  and  (2)  a three  nearest 
neighbor  decision  rule.  Each  texture  was  classified  as  a point  in  6-dimen- 
sional  feature  space.  For  the  first  measure,  a mean  and  standard  deviation 
in  each  dimension  were  calculated  for  each  texture  from  the  training  samples 
and  the  distance  for  an  unknown  sample  from  the  test  set  was  measured  from 
the  mean  in  standard  deviation  units.  The  results  are  shown  in  Table  2.  Also 
included  in  the  table  are  results  using  another  technique  discussed  in  the 
next  section. 

The  classification  matrix  for  the  method  using  the  3~nearest  neighbors 
decision  rule  is  shown  In  Table  3.  The  most  common  confusions  using  the  max- 
min  method  occur  among  wood,  fur,  and  water,  and  between  paper  and  cork.  It 
seems  that  the  distance  measure  created  by  the  algorithm  Is  similar  to  that 
used  by  humans  in  grouping  textures. 
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III.  COMPARISON  WITH  SPATIAL-DEPENDENCE  TECHNIQUE 


The  most  common  texture  classification  techniques  use  statistical  mea- 
sures based  on  spatial  dependence  probabilities.  In  order  to  make  the  tech- 
niques comparable,  we  limited  the  technique  described  by  Harallck  et.al.  [2] 
to  6 features  and  to  one  dimension.  We  used  64  grey  levels  in  the  spatial 
dependence  matrix.  The  best  results  we  could  obtain  in  using  this  technique 
was  83^  accurate  on  training  samples  and  66%  accurate  on  test  samples  as  shown 
in  Table  2.  The  classification  matrix  for  the  result  is  shown  in  Table  4. 

With  14  features  in  each  of  two  dimensions  the  spatial  dependence  technique 
accuracy  was  increased  to  94%  accurate  on  training  samples  and  89%  accurate  on 
test  samples.  However  the  computation  times  on  the  CDC  6500  were  indicative 
of  the  relative  efficiency  of  the  Max-Min  Technique.  The  total  feature  gener- 
ation and  classification  for  49  samples  of  each  of  8 textures  required  250  cpu 
seconds  for  the  6 feature  Max-Min  Technique,  845  cpu  seconds  for  the  6 feature 
spatial  dependence  technique,  and  1690  cpu  seconds  for  the  28  feature  spatial 
dependence  technique.  However,  it  should  be  noted  that  little  effort  was  made 
to  minimize  the  running  time  of  any  of  the  above  techniques. 

IV.  DISCUSSION 

The  max-min  appears  to  be  a promising  measure  of  texture  characteristics. 
The  method  has  presently  been  used  in  one  dimension  only.  One  two  dimensional 
extension  would  be  to  measure  the  max-min  features  in  several  directions  and 
use  that  direction  which  maximizes  some  criteria  plus  the  orthogonal  direction. 
This  would  make  the  algorithm  invariant  to  texture  rotation  as  well  as  adding 
a measure  of  the  rotational  symmetry  of  the  textures.  The  method  might  also 
be  used  to  detect  texture  boundaries  by  measuring  features  In  two  opposite 
directions  from  a suspected  boundary. 
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It  is  fairly  easy  to  incorporate  the  computing  required  for  the  max-min 
feature  extraction  in  special-purpose  hardware.  This  would  make  real  time 
texture  analysis  possible.  This  is  very  important  for  applications  such  as 
steel  mill  output  monitoring  where  a decision  must  be  reached  quickly  as  to 
whether  to  let  the  metel  continue  cool i ng  or  to  reprocess  lt» 

Also  the  quantities  measured  here  (number  of  extrema  vs.  threshold)  might 
be  called  a first  order  effect.  The  two  curves  for  cork  and  paper  in  Fig.  6 
are  almost  identical  and  confusions  might  necessarily  be  expected  in  a clas- 
sification algorithm  which  uses  only  slopes  of  these  curves.  However,  second 
order  measurements  which  include  information  as  to  how  the  small  extrema  are 
interspersed  among  the  large  extrema  would  differentiate  between  these  two 
textures . This  might  be  the  beginning  of  a hierarchical  structure  of  texture 
primitives:  those  that  differ  in  first  order  measurements  and  those  that 

differ  in  second  order  measurements. 
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Table  1 Training  Set  Data  Showing  the  Feature 
Values  for  a Few  Texture  Samples 


Features 


Name  of 
Texture 

1 

#@130 

2 

ii  @ 110 

3 

# @ 90 

4 

# @ 70  ’ 

5 

a e 50 

6 

# e 30 
a @ 10 

a @ no 

4 e 96 

a % 70 

# e 50 

a @ 30 

D4 

Cork 

.7621* 

.7637 

.7338 

.5077 

.6246 

.6793 

D4 

Cork 

.7126 

.7767 

.6782 

.6788 

.6227 

.6287 

D4 

Cork 

.*♦957 

.6964 

.7671 

7821 

.7467 

.5327 

D70 

Wood 

0 

.0385 

.2653 

.3793 

.5228 

.6520 

D70 

Wood 

0 

.0405 

.3318 

.3940 

.5345 

.7209 

D70 

Wood 

0 

0 

.0641 

.3805 

.3721 

.4398 

D57 

Paper 

.5943 

.7465 

.7513 

.8008 

.7173 

.4301 

D57 

Paper 

.6182 

.6875 

.7619 

.6383 

.5557 

.5267 

D57 

Paper 

.5319 

.6812 

.7886 

.7479 

.7112 

.4202 

D9 

Grass 

.6373 

.7751 

.8328 

.7456 

.7817 

.6636 

D9 

Grass 

.6854 

.6893 

.7536 

.7765 

,7447 

,6904 

36 
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Table  2 Classification  Results  Using  36 
Training  Samples  and  13  Test  SarT 
Each  of  8 Texture  Patterns 


Features  Used 


1.  Max-Mln  Method 

Extrema  Ratios,  6 Features 


2.  Max-Min  Method 

Normalized  Extrema  Ratios, 
6 Features 


3“Neai  ist  Neighbor 
Decision  Rule 

Training 


Weighted  Distance 
Decision  Rule 


Training  Test 


3.  Spatial  Dependence  Method 
6 Features,  One  Dimension 


Spatial  Dependence  Method 
28  Features,  Two  Dimensions 


Table  k Classification  Matrix  for  Spatial  Dependence  Method 

Using  6 Features  and  3-Nearest  Neighbor  Decision  Rule 


Texture 


Di* 

Cork 

1 

D70 

Wood 

2 

D69 

Wood 

3 

D93 

Fur 

k 

D29 

Sand 

5 

D57 

• 

Pap6r 

6 

D38 

Water 

7 

D9 

Grass 

8 

Assigned  Category 

Training  Samples  Test  Samples 

1 2 3 4 5 6 7 8 | 1 2 3 4 5 6 


---23 
3 - - - 
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APPLICATION  OF  OUTER-PRODUCT  EXPANSIONS  TO 
FEATURE  EXTRACTION  IN  PICTURES 

K.  Fukunaga  and  G.  V.  Sherman 

A brief  description  of  the  two-dimensional  version  of  the  subspace  method 
follows.  For  simplicity  square  nxn  pictures  will  be  assumed.  First  a 
representation  space  is  constructed  for  each  of  M classes.  A typical  repre- 
sentation space  is  characterized  by  the  orthogonal  projection  operation  Q(») 


Q (A) 


i=l  j=l 


(4>T  a j ) r 


(i) 


where  A is  a random  nxn  picture,  and  m ^ n,  m < n.  m and  m are  the  column 

c r — c r 

space  and  row  space  dimensionalities  respectively  and  are  selected  so  that 
about  90S  to  95%  of  each  column  or  row  squared  norm  is  contained  in  the  repre- 
sentation space  on  the  average.  The  { } and  } 1=1 ,n  are  arbitrary 

orthonormal  bases.  However,  two  bases  in  particular  have  been  found  to  work 
well.  The  first  Is  the  two-dimensional  Karhunen-Loeve  transformation  [l]  in 
which  the  { <f>  .}  are  eigenvectors  of  the  sample  autocorrelation  matrix  obtained 
from  all  columns  of  training  pictures  In  a particular  class.  Similarly  the 

are  derived  from  the  sample  autocorrelation  of  picture  rows  from  the 
same  class. 

The  other  basis  referred  to  is  obtained  by  selecting  the  { <f>.  } and  { 'p  i 

i i 

as  the  singular  vectors  of  the  singular  value  decomposition  of  E { A}  [2]  where 
expectation  E is  restricted  to  the  class  In  question. 

Feature  spaces  Q'  are  obtained  from  the  representation  spaces  Q by  either 
(2)  or  (3).  Subscripts  Indicate  class  labels. 


This  work  was  supported  in  part  by  NSF  Grant  GJ-35722. 
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The  intersection  of  two  subspaces  is  defined  to  be  the  subspace  contained  in 
both  of  the  original  two  subspaces  [3].  The  union  of  two  subspaces  is  defined 
to  be  the  subspace  consisting  of  all  linear  combinations  of  vectors  original- 
ting  from  the  original  two  subspaces.  Note  that  we  are  dealing  with  the 
vector  space  of  nxn  arrays. 

Our  research  has  shown  that  (2)  and  (3)  are  definitely  not  equivalent. 
Part  of  this  phenomenon  is  due  to  failure  of  subspaces  to  satisfy  the  dis- 
tribution law  of  logic.  Namely 

Qi  n (<*2  u Q3)  4 (Q,  n o 2)  u (q1  n q ) (1*) 


u (QjH  Q3)  4 (Q1  U 02)  n (Q1  U Q3)  (5) 

However,  further  research  is  needed  to  study  this  phenomenon.  Algorithms  for 
calculating  U and  Qj  D Q2  exist  but  will  not  be  discussed  for  brevity. 
Classification  is  accomplished  by  minimizing  criterion  /'(•). 


m m 
c r 


J*(k  ) = min  E E *(k)  (A-£  { A} ) * fk) 

k i=l  j=l  1 J 


-*■  A £ U), 


If  class  minimizes  J (•)  over  all  other  classes  to^  then  random  picture  A 

is  classified  into  class  to  . 

ko 

Experiments  with  this  method  on  the  24x24  Munson  handwritten  numerals 
demonstrate  a phenomenal  feature  extraction  capability.  A 75%  reduction  in 
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the  "”bcr  °f  futures,  from  576  to  166,  results  in  a mere  2. 5*  decrease  In 
classification  accuracy. 

Further  research  is  needed  to  relate  the  intersection  of  feature  sub- 
spaces to  intrinsic  dimensionality,  data  dimensionality,  number  of  classes, 
and  classification  error.  The  power  of  this  method  lies  in  its  reliance  on 

data  structure  determined  by  training  samples  rather  than  dubious  parametric 
forms . 

Preliminary  investigations  in  these,  areas  are  promising.  For  instance 
we  already  found  that  M subspaces  always  intersect  if  the  sum  of  their 
dimensionalities  is  strictly  greater  than  (M- 1 ) times  the  dimensionality  of 
the  observation  space.  So  the  problems  are  solvable. 
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A BRANCH  AND  BOUND  CLUSTERING  ALGORITHM 
K.  Fukunaga  and  P.M.  Narendra 

In  clustering,  each  of  N samples  is  assigned  to  one  of  clusters  w j , . . . 
where  the  number  of  clusters,  M,  is  assumed  to  be  prespecified.  Then,  there 
are  MN  possible  assignments  and  they  are  expressed  in  a tree  form  as  follows: 

SAMPLE 

X, 

z 

*3 

1 

Fig.  I 

In  order  to  apply  the  branch  and  bound  method  successfully,  we  have  to  come  up 

with  a proper  criterion  to  evaluate  each  node  so  that  the  branches  under  the 

node  can  be  eliminated  from  the  search. 

One  of  the  popular  criteria  which  has  been  used  in  clustering  frequently 

is  J = tr  S * S where  S and  S are  the  mixture  scatter  matrix  and  the  within 
m w m w 

class  scatter  matrix  respectively.  Since  the  coordinate  can  be  selected  to 

satisfy  S = 1 , we  may  use  J = tr  S . This  criterion  was  first  used  in 

ISODATA  [1],  [2].  Thus,  our  problem  is  to  find  the  cluster  assignments, 

c. so  as  to  minimize  tr  S . 

I N w 


ASSIGNMENT 
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When  the  search  comes  to  a node  at  the  k-th  level,  C) ck  are  already 

given  (for  example,  Cj  = and  c^  = for  the  node  A of  Fig.  1),  but 

Ck+1 ***** cn  cou,d  be  anything.  Therefore,  if  we  can  calculate  either  ck+1>?l?,c 

J(C|,.,.,ck,  ck+l CN^01*  ^ower  bound  of  that 


L(C| ,..., ck ) 


<__  min 
Ck+1 


J (c1 .... ,ck,ck+) , ... „cN) 


(1) 


then  with  the  satisfaction  of 


B 


L(e. 


(2) 


all  branches  under  the  node  can  be  eliminated  from  the  search  where  B is  the 
current  lowest  J found  up  to  the  present.  That  is,  all  possible  cluster  assign 
ments  under  the  node  give  larger  J's  than  the  one  which  was  already  found. 

The  tighter  the  lower  bound  Is,  the  more  branches  are  eliminated  effectively 
from  the  search. 


Our  study  levealed  that  the  wlthln-class  scatter  criterion  has  a nice 
property  as 

N J(cl CN)  -k  J(c! ck>  + (N'k)  J(ck+1 V 

where  J(cJt...,ck)  is  calculated  from  the  within-class  scatter  matrix  of 

X1 \ with  c,uster  assignment  c,,...,^,  and  J (ck+J cN)  is  for 

V * 

k+1 aN  Wlth  ck+l CN*  Let  J ^ck+l CN^  be  the  criterion  vali  for 

the  optimum  cluster  assignment  when  only  a subset'  { XR+] XN)  is  cc.^idered. 

Then,  J (^k+j c^)  <_  J (ck+j , . . . ,c^) , and  (3)  Is  bounded  from  the  lower 

side  as 


N J(c),...,cN)  >_  k J(Cj ck)  + (N-k)  J*(ck+J cN)  (A) 

The  right-hand  side  of  (k)  Is  Independent  of  the  cluster  assignments  for 
^k+l,****^N»  and  provides  the  lower  bound  of  (1). 


P * 

% 


! 


There  are  many  possible  ways  to  calculate  J (ck+? , . . . ,cN) . One  is  to 
start  from  J (cN)  and  to  calculate  J (c^,^),  J (CN.2»  cn-1‘  success- 

ively. The  branch  and  bound  method  with  the  node  evaluation  of  (4)  gives  an 
efficient  way  to  do  this  successive  extension.  When  the  number  of  samples 


becomes  large,  we  may  divide  the  samples  into  L groups  ,'{  X ,...  ,X,  }, 

1 k] 

* *k  ^ +j,...,Xk  ^ (^|_  = N)*  The  cluster  assignment  of 

• 2 L- 1 L 

each  subset  is  optimized  as  was  mentioned  above,  keeping  the  results  of  the 


. ..  h-v-v  ‘•'■'■iuiu,  « /,  o /,...  . men,  two  suDseis, 

i i i 

for  example  {Xj,,..,Xk  } and  { X^  +j,...,Xk  },  are  combined  to  forma  new  la 

• 2 

larger  subset.  For  1 <k  <k^,  the  lower  bound  is  computed  by 


(t<i  + k^)  J(c]  ,...,ck  +k  ) >.  k J(cj,.,.,ck)  + (kj  k)  J (ck+|,...,ck  ) 

+ (k2  - k j ) J*(c^+1 c(  ) (5) 

Thus,  L subsets  are  grouped  to  larger  L/2  subsets,  and  they  are  again  grouped 
to  LA  subsets  and  so  on. 

The  procedure  to  calculate  the  lower  bound  mentioned  above  is  only  one 
example  among  many  possibilities.  We  would  like  to  try  many  others.  Also,  we 
would  like  to  apply  this  idea  to  other  criteria,  particularly  to  the  valley- 
seeking clustering. 
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image  restoration  using  the  projection  method  algorithm 

T.S.  Huang,  M.  Kaveh  and  S.  Berger 

' • Introduction 

pr°Jecti°"  has  shown  promise  as  an  effective  method  of 

'ma9e  reSt°rat'°'’-  l(  ,s  -I, -suited  to  the  treatment  of  ilnear  spatial, y- 
variant  degradations.  A, so,  certain  types  of  a priori  Information  that  may  he 

aval iable  about  the  original  image  can  be  easliy  incorporated  into  the  res, or- 
ation  process. 

The  algorithm  has  been  applied  to  both  one  and  two-dimensional  signals  in 
order  to  evaluate  Its  potential  as  a restoration  technlgue.  The  results  which 
are  outlined  In  this  report  indicate  that  the  algorithm  can  be  an  effective 
method  for  the  restoration  of  images. 

2.  The  Algor i thm 

If  the  degradation  of  a two-dimensional  Image  can  be  represented  by 

9(*,y)  =■  D[f(x,y)J  +n(x,y),  where  f(x,y)  Is  the  original  Image,  n(x,y)  is 

noise.  0 is  a degrading  operator,  and  g(x.y)  Is  the  degraded  image,  then  the 

purpose  of  the  restoration  process  Is  to  approximate  the  original  Image  as 

closely  as  possible.  The  degrading  operator  is  assumed  to  be  linear  although 
space-variant. 

If  we  neglect  noise,  the  discretized  version  of  the  degradation  is: 

91  = al  I f 1+a12f2+*  “+aiwf 


92  a2|fi+  •••  +a?wf 


IN’  N 
2NfN 


9M  ~ aMlfl+  •“ 


+aMNfN 


Where  N may  not  egual  H.  The  projection  algorithm  Is  an  Iterative  technlgue 
for  solving  this  system  of  equations. 


N 


The  solution  is  obtained  by  successive  projections  onto  planes  in  hyper 


space.  For  the  case  where  a unique  solution  does  exist,  the  algor i thm  wi 1 1 
converge  to  the  point  of  intersection  of  the  hyperplanes.  If  the  planes  do 
not  intersect  at  a single  point,  the  algorithm  will  converge  to  a point  which 
may  be  a useful  approximation  in  a restoration  sense. 

The  mathematical  form  of  the  algorithm  is  given  in  vector  form.  Let 
f_  = (f  j , f , . . . , f |^)  be  a point  in  N space.  If  the  initial  guess  at  a solution 

• f (o  ) r (o  ) r (o)  r \0 ) , « • • « • • r . i 

is  f. = f »^2  * • • • »' fj  » the  new  solution  is  the  projection  of  the 


point  f 


(o) 


onto  the  hyperplane  g^  = al 1 ^1 +* * * 43  1 N^N*  *'° 


,(D  . f(o)  £(°),v9i1  Sj 

— i *— i 

where  a^  = (a^,a^ a|N^*  nex*  iteration  consists  of  the  projection 

of  f_^  onto  the  hyperplane  = a^  • f_  where  ~ (a2 ^ , . . . ,a2^) . This  pro- 
cess is  continued  to  the  Mth  plane.  This  completes  one  cycle  of  iteration. 


3.  Results 

The  effectiveness  of  the  projection  method  has  been  evaluated  for 
several  test  casesl  One  test  involved  a spatially  variant  degradation  in  the 
horizontal  direction.  The  original  test  image  was  an  "X"  in  a 128x128 
pictule  array.  The  intensity  of  each  picture  element  was  represented  by  an 
integer  from  0 to  25T.  The  image  arrays  were  stored  on  magnetic  tape  in  a 
compacted  format.  The  Gould  electrostatic  plotter  was  used  to  obtain  half- 
tone reproductions  of  the  images.  The  restoration  was  implemented  on  a CDC 
6500  computer. 

The  degradation  was  a smearing  which  increased  its  effect  toward  the 
edges  of  the  image.  The  points  on  the  vertical  center  line  were  not  affected. 


By  applying  the  projection  algorithm,  the  Image  was  made  Increasing,,  sharper 

With  each  iteration.  However,  the  effect  of  any  noise  was  also  increased. 

So  the  user  would  have  to  decide  which  Iteration  yielded  the  "best" 
restoration. 

The  performance  of  the  projection  method  was  also  investigated  for  a one- 
dimens iona,  signal.  The  original  signal  consisted  of  two  unit  pu,Ses  separa- 
ted by  a small  distance.  The  degraded  signal  was  obtained  in  the  fourler 
transform  domain  by  the  use  of  a triangular  ,™-pass  filter.  The  resultant 
degradation  caused  the  two  pulses  to  be  smeared  into  one  wide  pulse. 

The  projection  algorithm  is  capable  of  utilising  a prior,  information 
ehout  the  s.gnal.  F„r  instance  if  the  signal  is  known  to  be  positive  or  be 
tero  certain  regions,  the  algorithm  can  Include  this  knowledge  In  the 
restoration  process.  Sev.-rla  cases  were  treated  with  various  amounts  of  a 
Priori  information.  The  aigorithmwas  capable  of  resolving  the  image  into 
distinct  puises  after  a few  iterations.  In  general,  the  addition  of  more 
e Priori  data  improved  the  performance.  The  projection  rethod  compared  favor- 
-ly  with  the  least-squares  inverse  filter  method,  which  was  also  employed. 
Two-dimensiona,  images  of  ,28x128  points  were  also  tested.  The  Images 

were  dig.tiaed  versions  of  photographs.  The  results  Indicate  that  the 
algorithm  can  improve  on  thP  • 

e on  the  degraded  image,  but  the  effects  of  noise  dominate 

the  image  after  a certain  number  of  iterations  f u 

Keratlons-  Each  complete  iteration  re- 
quires about  12  seconds  of  computation  time  Thf>  *• 

con  time.  The  opt  .mum  number  of  i terat  ions 

is  a subjective  quantity. 

All  the  above  experiments  Invoived  Image  degradations  which  were  simulated 
computer.  In  our  recent  work,  we  applied  the  projection  algorithm  to  an 
optically  degraded  Image.  Specially,  we  took  an  image  degraded  by  camera  motion 
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(Fig.  1(a)),  digitized  it,  and  fed  it  into  the  computer.  Then  it  is  restored 
by  Wiener  inverse  filtering  and  by  the  projection  algorithm.  The  results  are 
shown  in  Figs.  1(b)  and  (c).  The  projection  algorithm  gave  a much  better 
restored  imacje. 
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ESTIMATING  THE  IMPULSE  RESPONSE  OF  A DEGRADING  SYSTEM 


Brian  O'Connor  and  T.  S.  Huang 

One  major  problem  in  image  processing  is  the  estimating  of  the  impulse 
response  of  a degrading  system.  Once  this  is  found  a inverse  filter  or  other 
algorithm  can  be  found  and  applied  to  the  image  to  ci;  nate  the  effects  of  the 
degrading  system.  The  estimation  problem  car  be  broker,  down  into  two  cases; 
in  the  first,  the  ideal  and  degraded  images  are  assumed  to  be  known  to  some 
extent;  and  in  the  second,  only  the  degraded  picture  is  present.  This  latter 
case  is  called  blind  deconvolution.  Introduced  below  is  a technique  which  can 
be  applied  to  either  of  the  above  cases,  but  should  find  wider  use  in  the 
latter.  This  method  incorporates  image  segmentation  with  a modification  of 
Knox's  method  [ 1 j for  multiframe  processing  to  estimate  both  the  magnitude  and 
phase  of  the  degrading  system  function. 

Let  v(x,y)  be  an  ideal  image  which  Is  degraded  by  a system  with  impulse 


response  h(x,y).  The  observed  image  is  vc(x,y)  = v*h  which  can  be  segmented 
into  many  rectangular  strips,  vc.(x,y),  i = 1,...,N.  If  the  extent  of  h(x,y) 
is  small  compared  to  the  duration  of  the  segments  vc.(x,y),  then  vc.(x,y): 


v i " h • This  implies  VC.^,^)  = V.  (f  j ,f£) -H(f  ] J ; furthermore,  |VC.  (f)  ,f2) 


- |vr(f,.f2)|-|H(f1,f2)|  and  /_  VC.  (f,  ,f£)  -^(f,,^)  +^H(f),f2).  To 


estimate  | H (f  1 f^\  the  square  magnitudes  of  the  N segments  are  averaged,  so 


I H (f . , f „ 


^ TT  Elvc 


i 


In  some  applications  the  light  distribution  of  the  N regions  of  the  original 
image  vary  sufficiently  fast  and  are  sufficiently  different  from  one  another 

I 2 

so  that  N l |V(f,  , f 2 )|  is  approximately  a constant.  This  Implies  that  the 

approximation  of  the  magnitude  Is  Independent  of  original  scene. 
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order  to  estimate  the  phase  of  the  system  function  an  average  auto- 
correlation  is  performed. 

1 VCi (fl ’V'VCj (f|+Afl »f2+AV 


[£  vi  (f,  *f2)*,ri  (f,+*f, .f2+af2)]  H(f,  ,f2)*TT{fl+af)  ,^+af^) 
After  some  manipulations  we  find 


(f, .f2)-vc,  (f ,-faf | . f ?+af ?) 1 1 ry . (f | ,f;).y,  (f^af,  ,f„+af„)  | 

I ^c, < f t > f 2 ) * vc , (f,,4f,,f2+af2)|  [rv| (f f ,f2)-v. (f(+af ) ,f2+af2)] 
H(f..f2)-H(f  +af  f +1f  ) 

' ~ T , ~ — ~ - ej(e(fi-f2)  - e(fI+Af].fz^f2J 


l"(f,.f2)'H(f|+4f  f +4f  )| 


that  the  imaginary  part  of  the  comp, ex  logarithm  gives  an  estimate  of  the 
Phase  in  terms  of  phase  differences.  In  the  above  manner  the  phase  differences 
between  every  pair  of  adjacent  points  can  be  measured  over  the  whole  transform 


plane. 


Programs  are  being  written  to  simulate  the  above  technique  and  should  give 
better  results  than  procedures  which  estimate  only  the  magnitude  of  the  de- 
grading system  function.  A problem  could  arise  in  the  phase  estimation  since 
it  is  determined  by  add.ng  phase  differences  together,  which  means  that  the 
Phase  errors  add,  thus  producing  the  variance  of  the  final  error  to  increase 


1 inearly. 
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MISSILE  TRACKING  ALGORITHMS 


T.  S.  Huang  and  G.  Y.  Tang 

This  project  is  motivated  by  the  real-time  video  tracking  problems  as 
the  U.S.  Army  White  Sands  Missile  Range.  The  missile  movie  frames  were 
supplied  to  us  by  Dr.  Alton  Gilbert  of  WSMR. 

Our  preliminary  experimentation  has  indicated  that  it  is  possible  in  many 
cases  to  detect  and  locate  a missile  by  processing  single  scan  lines.  More 
specifically,  we  have  found  that  it  is  possible  to  determine  whether  a scan 
line  passes  through  a missile  (and  the  location  of  the  crossing)  by  cross- 
correlating  the  scan  line  with  a paradign.  This  opens  up  the  possibility  of 
very  fast  algorithms  which  scan  selectively  instead  of  full  frame.  Solid  state 
imagors,  such  as  surface  acoustical  wave  and  charged-transfer  devices,  are 
particularly  suitable  for  selective  scanning.  These  devices  can  also  be  used 
to  do  cr 3sscorrelation  at  extremely  high  speed. 

1 ) Missile  Tracki ng 

We  describe  an  example  of  missile  tracking  using  single  scan  line  cross- 
correlations. We  used  10  evenly  spaced  scan  lines  over  the  missile  frame  shown 
in  Fig.  1,  and  crosscorrelate  each  scan  line  with  the  paradign  shown  in  Fig.  2. 
Bv  applying  a suitable  criterion  to  the  correlation  results,  we  were  able  to 
determine  which  scan  lines  passed  through  the  missile  and  the  crossing  loca- 
tions. Then,  among  the  scan  lines  that  passed  through  the  missile,  we  picked 
out  the  uppermost  and  lowermost  scan  lines;  and  we  scan  at  a finer  spacing 
than  before  a few  lines  above  the  uppermost  and  a few  lines  below  the  lower- 
most to  locate  the  noise  and  the  tail  end  of  the  missile  (again  by  single  scan 
line  correlations).  The  resuU  is  shown  in  Figs.  3 and 

2)  Initial  Acquisition  of  the  Missile 

in  the  WSMR  application,  we  know  beforehand  where  the  missile  will  be 


'N 
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launched.  We 


doing 

site. 

1 ines , 


can  aim  our  imagor  at  the  ieunch  site  and  scan  repeatedly  (and 
crosscorrelation)  at  high  speed  one  or  several  lines  above  the  missile 

Then,  we  can  detect  the  missile  as  soon  as  It  moves  through  these  scan 
and  try  to  track  it  from  that  point  on. 


INFORMATION  EXTRACTION  FROM  y-RAYS  IMAGES* 

G.  Y.  Tang 

Gamma  ray  cameras  have  been  used  by  physicians  for  many  years  fo  detect 
diseases.  The  grey  tone  distribution  of  the  picture  obtained  by  the  gamma  ray 
camera  contains  most  of  the  information  about  the  mass  density  of  the  patient's 
'issues  or  muscles.  As  a matter  of  fact,  the  gamma  ray  picture  is  a projec- 
tion of  the  mass  density  of  the  object  illuminated  by  a radiation  source. 
Physicians  make  their  judgement  by  comparing  an  unknown  picture  with  some 
known  diseased  patterns.  However,  since  the  rich  varieties  of  diseased  pat- 
terns and  the  mutual  influence  between  symptoms,  the  procedure  of  making 
judgement  is  not  a straightforward  work.  Usually  it  takes  two  or  three  years 
training  plus  experience.  In  some  cases  it  Is  even  required  a joint  judge- 
ment of  several  physicians.  The  complexity  in  judgement  does  not  mean  that 
the  request  from  patients  should  be  ignored.  A computer-aided-judgement- 
making  algorithm  is  therefore  necessary  in  order  to  alleviate  the  physician's 

burden  and  to  serve  more  patients.  This  is  also  the  ultimate  goal  of  studying 
gamma  ray  pictures. 

m this  report,  an  attempt  has  been  made  to  employ  pattern  recognition 

techniques  to  solve  the  problem  of  compute r-a i deded  judgement  making  for  gamma 

ray  pictures.  As  shown  in  Fig.  I,  a pattern  recognition  system  consists  of 

three  stages  in  sequence,  i.e.,  pattern  analysis,  feature  selection  and 

classification.  Pattern  analysis  is  to  see  what  is  the  most  informative  part 

which  can  be  obtained.  Feature  selection  is  to  throw  away  some  unimportant 

measures  and  to  retain  only  those  which  are  sufficiently  representative  so 
that  the  computational  effort  in  classification  is  reduced  and  that  the 
probability  of  wrong  classification  due  to  the  disturbance  of  unrelevant  data 

<s  reduced  too.  The  purpose  of  classification  Is  to  assign  an  attribute  for 
each  input  pattern. 


This  work  was  supported  in  part  by  Consiglio  Nazionale  delle  Ricerche,  Italy. 
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Some  pattern  analysis  experiments  whcih  have  been  done  are  described  in 
the  first  section  of  this  report.  In  the  second  section  some  possible  ap- 
plications of  these  measurements  for  the  purpose  of  combined  aid  system 
recognition  are  discussed  briefly. 

I.  Pattern  analysis. 

A gamma  ray  picture  is  represented  by  a matr  ix  / (I  (i  , j) } i=I,...,N  , 

j*l,...,N  where  and  are  the  spatial  extention  of  the  gamma  ray  picture 

in  x and  y directions.  I ( i , j ) is  the  (i,j)  entry  of  that  matrix.  Its  value 

corresponds  to  the  intensity  of  radiation  energy  at  (i,j)  position.  In  our 

experiment,  N = N = 6**,  I ( i » j ) varies  between  0 and  2^-1.  Statistically  a 
x y 

gamma  ray  picture  is  an  outcome  of  a Poisson's  process.  So  P(l(i,j)  = x)  = 


■m  ( i , j ) 


m ( i ,j)' 


P(l(i,j)  = x)Ax  is  the  probability  of  I ( i , j ) of  value  x. 


x ! 

Poisson's  distribution  indicates  that,  for  the  area  of  higher  mean  value,  the 
variance  or  the  noise  is  larger  too.  In  order  to  combat  the  noise  disturbance, 
a smooth  technique  is  required.  A two  dimensional  symmetric  low-pass  filter 
is  cnoson  .'or  this  purpose.  After  the  noise  has  been  cleared  out,  we  threshold 
the  picture  by  a preset  threshold.  The  thresholded  picture  has  only  two  grey 
tones  i and  0.  1 is  assigned  to  those  picture  points  with  I ( i , j ) >_  T,  T is  the 

threshold,  0 is  otherwise.  The  thresholded  picture  is  denoted  by  l'(i,j),  A 
simple  edge  detector  can  be  applied  to  I " ( i , j ) in  order  to  get  contours.  The 
edge  detector  works  as  if  | !"(•,))  - l'(i+l,j)|  > 0 then  It  reports  an  edge 
point  at  (i,j),  otherwise  it  reports  no  edge  points.  Fig.  2 shows  the  contour 
obtained  by  setting  threshold  T=8000  and  by  using  the  foregoing  edge  detector. 
Notice  that  the  left  lobe  and  the  right  lobe  are  not  distinguishable.  This 
gives  u'  a hint  to  try  a iower  threshold  T=^500.  Fig.  3 displays  the  contour 
thus  obtained.  Both  lobes  are  not  shown.  A third  trial  is  a threshold  of 
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value  7000.  Pig.  k displays  the  result.  We  can  see  that  the  two  lobes  are 
present,  but  the  left  one  is  rather  small.  More  than  that,  the  deformity 
which  contains  most  of  the  symptomatic  information  is  not  significant.  All 
these  evidences  indicate  the  difficulty  to  set  up  a good  threshold  so  that  an 
informative  contour  may  be  obtained.  In  order  to  combat  that,  we  proposed  an 
approach  which  should  achieve  that  1)  it  guarantees  two  separate  contours  for 
the  two  lobes  2)  it  subjects  to  a better  tolerance  for  the  variation  of  the 
th reshol d (Th i s means,  for  example  T=^500  and  T=7000  should  generate  about  the 
same  contours).  The  basic  frame  for  the  proposed  method  is  to  break  up  the 
entire  picture  into  several  regions.  Then  we  locate  these  two  regions  which 
enclose  the  two  lobes.  New  threshold  on  the  original  picture  for  each  of 
these  two  regions  is  set  up  according  to  the  local  grey  level  distribution. 
Finally,  a contour  follower  is  employed  to  get  the  coordinates  of  the  points 

on  the  contour.  The  foregoing  edge  detector  is  used  for  the  purpose  of  dis- 
playing. 

More  specifically,  the  proposed  method  can  be  described  as: 

(1)  Picking  up  a threshold  T. 

(2)  An  operation  T,  is  defined  as 

h 

V'  .j)  = 1 if  I (i,j)  > T 
Th(i  *j)  = 0 if  I (i  ,j)  <T 

and  I - Th(|)  means  l'(i,j)  = Th ( i , j ) . So  \'  is  a matrix.  I ( i , j ) which  is 
the  picture  matrix  has  been  defined  before,. 

(3)  An  operation  E,  which  is  the  edge  detector,  is  defined  as: 

= 1 " I'O+l.j'l  > 0 for  all  i and  j 

E(i,j)  = 0 otherwise,  where  I'  Th  ( I ) and  I"  - E(T)  means 
1 (i » j ) * E (i , j ) for  all  i ,j  . 
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(*0  From  \",  we  obtain  two  points  TOP,  BOTTOM  as 


suc^  = 1 and  jj  Is  the  maximum  of 

all  j such  that  E(i  j)  * |. 

BOTTOM  = (12J2)  such  that  E(i2fj2)  = 1 and  j2  is  the  minimum 
of  all  j such  that  E ( i , j ) = 1„ 

(5)  A line  M is  defined  as 

M = <(U)|  J = j (jj+jj),  i = I Nx> 

(6)  A line  M^  is  defined  as 

Mj  = {('»j)  lj  = J 02+  j (j,+j2),  i = 1 N } 

(7)  An  operation  C on  each  horizontal  line,  i.e.  j = const,  is  defined  as 

1 N 

C 0 ‘ COnst)  = — EX  r (?,j)  xi 

j M 

where  N.  is  the  number  of  l's  on  the  line  j = const  on  picture  \\  This 
operation  locates  the  centroid  on  each  horizontal  line. 

(8)  For  lines  lying  between  M and  M, , we  can  apply  operation  C to  each  of 

them.  Then  we  obtain  a set  of  H -H+l  points  (x  y > , (x  y ) (.  u , 

11  z L M. -M+l , 

^-M+P*  The  first  component  is  obtained  by  operation  C,  the  second  is  the 

one  defining  the  horizontal  line. 

(9)  A line  V = AX  + b is  found  such  that  ' Z ’ (y  -ax  -b)2  is  minimum. 

i = l 

V is  corresponding  to  the  second  element  and  X Is  corresponding  to  the  first 
element  in  our  l(l,j),  l'(l,j),  l"(i,j)  notations. 

do)  A line  segment  R:Y  - axtb  where  I (J,UZ)  < V<_j,  is  thus  obtained. 
(II)  The  picture  is  now  segmented  into  four  components  denoted  by  A,  B,  C, 
o.  A is  the  background.  B is  what  is  below  the  line  M but  not  the  background. 
C is  what  is  to  the  left  of  R and  what  Is  above  M,  but  not  the  background.  0 
is  the  same  as  C except  to  the  right  of  R,  Fig.  5 shows  the  segmentation  from 


N 


Fig.  2.  Fig.  6 shows  the  segmentation  from  Fig.  3.  Fig.  7 shows  the  seg- 
mentation from  Fig.  b. 

(12)  From  now  on,  we  are  only  interested  in  area  C and  D.  For  each  of 

these  two  areas,  we  define  a histogram  from  which  a new  threshold  can  i.  e 

obtained.  The  common  maximum  and  minimum  in  areas  C,  and  D are  located. 

They  are  denoted  by  MAX  and  MIN  respectively.  The  histogram  H and  H is 

c d 

defined  as 


Hc0<)  = # {(i,j)|k*<  L0jj)~MIN  X I CELL  < k+l  ,- (i  , j ) EC  } and 

MAX-M I N 

Hd(k)  « if  {(i,j)|kj<  COCELL  <k+l  , (i  ,j)  E 0} 

MAX-M IN 

where  k - 1 ,2 , . . . , I CELL,  and  ft  denotes  the  number  of  elements.  The  new  thres- 

holds  T and  T,  for  the  area  C and  D respectively  are: 

k 

T = min  { k | l h (i  ) > N } 
i-1  “ C 

k 

Td  = min  { k | E h ( i ) >_  Nd  } 

Nc  = $ { (i ,j) | ( i » j ) EC}  X PER  1 

Nd  = * Ui,j)  | (i  ,j)  EC}  X PER  2 

PER  1 and  PER  2 are  two  preset  numbers  which  are  0 <PER  1 <1  and  0 < 

PER  2 < 1. 

(13)  Thresholding  and  edge  detecting  same  as  (2)  and  (3)  are  applied  to 
areas  C and  D with  new  thresholds  T^  and  Trf  respectively.  (PER  1 = PER  2 = 0.^). 

Fig.  8 shows  the  contours  obtained  from  Fig.  2.  Fig.  9 shows  the  contours 
obtained  from  Fig.  3.  Fig.  10  shows  the  contours  obtained  from  Fig.  b. 

Notice  that  we  can  always  obtain  two  separate  contours  by  the  three 
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different  thresholds,  i.e.,  T = **500,  7000,  8000.  Also,  the  contours  from 
T = **500  and  T * 7000  are  exactly  the  same  and  the  contour  from  T = 8000  is 
very  similar  to  that  from  T = **500  and  T = 7000.  Especially  for  the  good 
lobes  [the  right  one],  they  are  almost  indistinguishable. 

(?**)  A contour  follower  is  designed  to  replace  the  edge  detector  so  that 
it  will  always  report  a closed  contour  no  matter  what  the  shape  of  the  contour. 
The  input  to  the  cor  tour  fol lower  is  a window  of  size  3x3.  The  center  of  the 
window  [0  on  Fig.  11]  is  already  known  on  the  contour.  One  of  the  eight 
neighbors  of  the  center  is  known  on  the  contour  too.  The  function  of  the 
contour  follower  is  to  find  a point  which  is  supposed  to  be  the  next  point 
in  the  eight  neighbors  around  the  center.  Then  the  window  moves  to  and 
centers  at  that  point.  Repeating  the  same  procedure,  we  car  obtain  tne 
whole  contour.  The  way  to  locate  the  next  point  from  the  eight  neighbors  of 
the  center  is  simply  looking  at  the  eight  neighbors  i t . counterclockwise  sense 
starting  with  the  neighbor  already  known  on  the  contour.  For  example,  on 
Fig.  IT,  searching  starts  with  2 and  then  3,**, 5 until  6 where  the  next  point 
is  found.  The  contour  then  reports  the  coordinate  of  the  next  point  (i,j) 
and  the  direction  of  moving  which  is  cod^d  into  8 numbers  as  shown  on  Fig.  11. 

Table  1 lists  the  two  contours  obtai  ted  by  setting  T = 8000  and  PER  1 = 

0.*4  and  PER  2 = 0.**.  The  edge  detection  is  replaced  by  the  contour  follower. 
The  corresponding  direction  code  is  attached  too.  Fig.  l*i  shows  the  scatter 
distribution  of  the  direction  code  on  Table  1. 

(15)  The  contours  obtainec  from  the  unfiltered  picture  by  the  same  method 
as  to  obtain  Fig.  10  is  shown  on  Fig.  13  where  T = 8000,  PER  1 = PER  2 = 0.**. 
Obviously  these  two  contours  are  more  noisy. 

Fig.  15  shows  the  contours  by  settiny  T « 8000,  PER  1 = 0.6,  PER  2 ■ 0.** 
from  a filtered  picture  by  the  same  method  as  to  obtain  Fig.  10. 
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Fiq.  16  shows  the  contour  by  setting  T = 8000,  P<ER  1 = 0.15,  PER  2 = Q.k 
from  a filtered  picture  by  the  same  method  as  to  obtain  Fig.  10. 

II.  D i scuss i ons . 

In  the  previous  section,  it  shows  the  possibility  to  obtain  two  separate 
contours  from  a gamma  ray  picture  by  the  proposed  method  with  which  a set  of 
coordinates  (x.,y.),  i=l,...,N  is  reported.  How  to  use  this  set  of  coordinates 
to  make  a classifier  is  another  problem.  One  way  to  tackle  this  problem  is  to 
perform  a t ransforma t i on  on  the  coordinates  so  that,  after  transformation,  the 
new  coordinates  should  be  1)  invariant  for  translation  of  the  object,  2) 
invariant  for  the  rotation  of  the  object,  3)  invariant  for  the  reflection, 

*0  invariant  for  the  swelling  or  shrinking  of  the  object.  Several  methods  ha 
have  been  studied  in  order  to  achieve  them.  Locating  the  centroid  and  using 
the  centroid  as  the  origin  of  a polar  coordinate  system  is  one  way.  The  use 
of  Fourier  trans forma t i on  by  treating  each  point  (X.,Y.)  as  a complex  number 

VJ'V  J'  = ,/-1  !s  another  way.  lhe  curvature  may  be  a useful  quantity  too. 
However,  whatever  the  transformation  is,  we  will  have  a new  set  of  coordinates 

i=1 M*  Generally  M is  smaller  than  N.  A discriminant  function 

f is  defined  on  u.,v.),  i = l M such  that  f (u]  ,v,  ,u2  ,v2 u,vV  - 0 

i nd i cate  disease,  and  f (u) , v| , u£ , v£ , . . . ,uM,uM)  >_  indicates  no  disease  f may  be 
a linear  or  a nonlinear  function.  Supervised  or  non-supervised  (learning 
techniques  may  be  used  to  get  f wi th  a set  of  training  data  which  we  know  if 
any  of  them  is  good  or  diseased. 

Accompanying  with  experiments,  further  studies  of  transformations  and 
discriminant  functions  are  necessary  to  complete  the  system.  The  syntactic 
pattern  recognition  technique  may  be  adapted  too.  How  acceptable  of  the 
syntactic  approach  depends  mainly  on  experimental  results.  It  is  hard  to 
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Fig.  2 Edges  obtained  by  setting  threshold 
equal  80CC. 


65 


I 

t 

I 


it  h ;*  » n » 

* HU  V*  N ft 


* N N 
III 


It 
I « 


I I 
ft  ft 


♦ ft 
«« 
ft 

«« 

ft 

ft 

ft 


ft 

ft 

ft  » 


ft  ft 
♦ ft 
ft  ft 
ft  ft 
ft  ft# 

* ft 
ft  ft  ft 
ft  ft  ft 

ft  ft  ft  ft  ft 


»t«*  tflti  until 


Fig.  3 Edges  obtained  by  setting  threshold 
equal  1»500. 
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Fig.  k Edges  obtained  by  setting  threshold 
equal  7000. 
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fig.  5 Four  regions  obtained  from  Fig.  2. 
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Fig.  7 Four  regions  obtained  from  Fig.  k. 
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Fig.  8 Contours  obtained  from  Fig.  2. 
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Fig.  9 Contours  obtaineu  from  Fig.  3. 
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Fig.  15  P contours  obtained  with  T=G000 
PER  1 = 0.6,  PER  2 = 0.4. 
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Fig.  16  Contours  obtained  with  T=8000, 
PERI  = 0.15,  PER  2 = 0.**. 
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SATELLITE  IMAGERY  NOISE  REMOVAL 
O.R,  Mitchell  and  P.L.  Chen 

I.  Introduction. 

This  recent  work  is  a continuation  of  the  previous  report  "Satellite 
Imagery  Noise  Removal"  [1].  We  are  using  homomorphic  filtering  techniques  to 
remove  multiplicative  noise  effects  such  as  cloud  and  atmospheric  turbulence 
in  ERTi.  imagery. 

Our  approach  is  to  estimate  the  noise  power  spectrum  by  classifying  each 
region  of  the  noisy  picture  according  to  the  level  of  the  noise  present  using 
the  mu  1 1 i spectra  1 data  analysis  software  system  developed  by  LARS.  Once  the 
noise  power  spectrum  is  obtained,  an  optimum  filter  is  derived  to  separate  the 
signal  and  noise.  Information  is  then  extracted  from  the  recovered  signal 
using  mul ti spectral  classification  techniques. 

II.  Theory  of  Cloud  Removal. 

Assume  an  image  of  the  earth  is  produced  when  a light  cloud  cover  exists 
over  the  region  of  interest  as  in  Fig.  1. 


Sun  Satellite 

Illumination  Scanner 
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We  have  assumed  that  the  cloud  reflection  of  sunlight  plus  the  cloud  absorption 
equals  one  and  that  the  illumination  is  approximately  constant  on  the  earth 
surface.  The  image  received  at  the  scanner  is 

s(x,y)  = alr(x,y)  t(x,y)  + l[l-t(x,y)]  < I 

where  r(x,y)  is  the  signal  and  t(x,y)  is  the  noise. 

sfx.y)  = It  (x,y)  [ar (x ,y)  -1]  + | 

Takino  the  logarith 

1 °9 [ |-s (x ,y ) ] - log  I + log  t(x,y)  + log [ l-ar (x,y) ] 

Not.ce  that  now  the  signal  and  noise  are  additive  so  that  an  optimum  linear 
filter  can  be  appl ied. 

The  technique  is  as  follows:  may  be  approximated  as  the  brightest  point 

in  the  image  [corresponds  to  t=0].  The  original  data  is  therefore  inverted  by 
subtracting  the  intensity  at  each  point  from  the  maximum  intensity  in  sfx.y). 

The  logarithm  is  then  taken  of  the  inverted  data.  Now  the  signal  and  noise 
are  additive.  The  power  spectrum  of  the  noise  t(x,y)  is  estimated  using  tech- 
niques discussed  in  the  next  section  and  a Wiener  filter  is  used  to  get  the 
estimate  of  1 og [ 1 -ar (x ,y ) ] . The  signal  is  then  recovered  by  taking  the  anti- 
log, inventing  the  intensities,  and  increasing  the  contrast. 

Ml.  The  Classification  of  Cloud. 

The  analysis  of  clouds  as  a multiplicative  noise  in  ERTS  imagery  is  im- 
plemented by  using  the  classification  technique  at  LARS.  In  this  report  the 

classification  of  cloud  is  based  on  the  amount  of  cloud  cover  located  at  each 
pixel. 

The  LARS  classifier  processes  mu! t i spectra  I data  one  point  at  a time 

classifying  unknown  data  using  training  statistics  developed  from  pre-classi- 
fied  data. 
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The  basic  assumptions  for  classification  are  as  follows: 

0)  All  class  probabilities  are  equal. 

(2)  The  probability  density  function  of  training  samples  is  Gaussian. 

Hie  decision  criterion  for  classification  is  the  so-called  maximum  like! 
hood  .ule.  The  unknown  data  point  is  compared  with  all  of  the  training 
classes  and  then  is  assigned  to  the  most  likely  class.  The  detailed  informa- 
tion about  classification  algorithm  can  be  found  in  LARS  information  notes 
[2,3]. 


A general  procedure  for  the  classification  of  clouds  is  described  in  the 
block  diagram  in  Fig.  2. 

Consider  Fig.  3.  This  is  a section  of  original  ERTS-1  mul t i spect ra 1 
scanner  data  (channel  A).  Obviously,  a large  number  of  classes  of  clouds  can 
be  defined,  however,  only  four  kinds  are  chosen.  These  are  Full  cloud  (F), 
Host  cloud  (M),  Half  cloud  (H),  and  Small  cloud  (.).  Water  and  ground  are 
defined  to  be  a class  of  no  cloud. 

Fig.  A shows  a classified  version  of  Fig.  3.  It  is  hard  to  compare  these 
two  pictures,  because  the  classified  picture  is  based  on  all  four  channels. 

The  description  of  cloud  as  we  defined  in  this  manner  has  proved  to  be 
very  successful  when  one  sees  the  training  and  test  field  performance. 

Tables  2-1,  2-2  show  these  results. 

IV.  Noise  Extraction. 

Fig.  5 shows  a small  section  of  Fig.  *4,  which  has  size  of  6*4x6*4  (line 
838-96*4-2 , Col.  1220-1  3*46-2). 

Notice  that  a classification  map  (Fig.  5)  ]s  simply  different  integers 
representing  each  class.  Since  these  are  six  spectrally  distinct  classes, 
irtegers  1 through  6 are  used  to  represent  them.  Fig.  6 shows  this  integer 
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N 


Group 


Table  2-1  Training  Field  Performance 


Number  of  Samples  Classified  Into 
No.  of  Pet.  Fulcld  Hoscld  Hahcld  Smlcld  Water  Ground 
Samps  Corct 


Fulcld 
Mosel d 
Hafcl d 
Smlcld 
Water 
Ground 


100.0  3 

100.0  C 
100.0  0 
100.0  0 
100.0  0 
100.0  0 


Overa! 1 
Performance 


G roup 


No.  of 
Samps 


38/  38)  = 100.0 

(From  LARS  *PR I NTRESULTS) 

Table  2-2  Test  Field  Percentage  (PCT.) 

Percentage  of  s 
Fulcld  Mosel d Hafcld  Smlr 


Percentage  of  samples  classified  into 
Hafcld  Smlcld  Water  Ground 


Fulcld 

10 

100.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Moscld 

2 

100.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Hafcld 

10 

0.0 

0.0 

80.0 

0.0 

0.0 

20.0 

Smlcld 

36 

0.0 

0.0 

0.0 

100.0 

0.0 

C.3 

Water 

8 

0.0 

0.0 

0.0 

0.0 

100.0 

0.0 

Ground 

6 

0.0 

0.0 

0.0 

0.0 

0.0 

100.0 

72 

16.7 

0.0 

11.1 

50.0 

11.1 

11.1 

Overa  1 1 

Performance 

68/ 

72)  = 3k. k 

(From  LARS  *PR I NTRESULTS) 
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F|g.  5 Classification  map  of  size  6*4x6*4. 

(-Lars  Run  73003600,  Line  13*46-1*47? 
Col.  1224-1350-2) . ' 
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Fig.  6 Integer  coded  classification  map 
of  Fig.  b. 
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classified  map.  A reasonable  guess  about  the  percentage 
fiected  energy  passing  through  the  cloud  is  listed  in  the 


transmission  of  re- 
fol lowing  table. 


Table  3-1  Cloud  blocking  information. 


Type  of  Cloud 

Percentage  of  Transmission 

Full  Cloud 

103' 

Most  Cloud 

30% 

Half  Cloud 

50% 

Small  Cloud 

75% 

Water 

100% 

Ground 

100% 

Consider  Fig.  6 and  Table  3“1  together,  a 
constructed  simply  by  transforming  each  coded  i 
percentage  value.  Fig.  7 shows  the  picture  for 


pattern  for  noise  only  can  be 
nteger  into  its  corresponding 
the  noise  only. 


fig.  7 Picture  of  cloud  (Half  and  Small)  and 
no  cloud  (Water  and  Ground). 


classified  map.  A reasonable  guess  about  the  percentage  transmission  of  re 
Mooted  energy  passing  through  the  cloud  is  listed  In  the  following  table. 


Table  3“1  Cloud  blocking  information. 


Type  of  Cloud 

Percentage  of  Transmission 

Full  Cl oud 

10% 

Most  Cloud 

30% 

Half  Cloud 

50% 

Small  Cloud 

75% 

Water 

100% 

Ground 

100% 

Consider  Fig.  6 and  Table  3"1  together,  a pattern  for  noise  only  can  be 
constructed  simply  by  transforming  each  coded  integer  into  its  corresponding 
percentage  value.  Fig.  7 shows  the  picture  for  the  noise  only. 


Fig.  7 Picture  of  cloud  (Half  and  Small)  and 
no  cloud  (Water  and  Ground). 
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V.  Homomorphic  Filtering  Process. 


The  two  dimensional  filtering  to  reduce  the  cloud  effect  is  accomplished 

by  weighting  each  point  in  the  discrete  Fourier  Transform  of  the  logarithm  of 

S. 

the  picture  grey  levels  by  a non  casual  filter  function  of  the  form  - ■ 

i 1 

where  S.  and  N.  are  the  vlaues  of  the  signal  and  noise  power  spectrums  at  the 
particular  frequency  point  being  considered.  Following  the  filtering  process 
we  take  the  inverse  DFT  and  then  exponentiate  to  obtain  the  filtered  picture. 
See  Fig.  8. 

VI.  Results  and  Conclusions. 

The  filter  function  tends  to  smooth  the  picture  resulting  in  a loss  of 
sharp  edges.  However,  the  synthesis  of  high-passed  picture  to  filtered  pic- 
ture will  compensate  this  disadvantage. 

The  overall  performance  of  the  filtering  processes,  described  in  this 
report,  will  be  evaluated  by  classifying  the  filtered  picture  and  comparing 
with  the  ground  true  map  of  the  same  area.  This  will  be  done  in  the  near 
future.  We  are  also  converting  the  filtering  operation  to  a three  dimensional 
filter,  to  take  advantage  of  the  four  spectral  channels  In  the  filter. 
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PDP-11  SOFTWARE  DEVELOPMENTS 


Janies  J.  Besemer 

The  major  accomplishment  on  the  PDP-11  this  year  was  the  installation  of 
a new  operating  system  called  UNIX.  It  was  developed  at  Bell  Labs,  and  is 
distributed  by  Western  Electric.  .he  University  of  Illinois  engineered  the 
system  modifications  to  allow  it  to  communicate  with  the  ARPANET.  Our  version 
of  the  system  has  been  fully  operational  since  early  October,  and  has  been 
operating  on  a 2k  hour  per  day  schedule  since  the  first  of  the  year.  We  are 
quite  pleased  with  the  new  system.  With  the  availability  of  UNIX,  most  of  the 
graduate  students  and  other  researchers  have  been  working  on  transporting  their 
software  from  other  sites  to  our  PDP-11.  Other  accomplishments  include  the 
(earlier)  installation  of  an  ELF  Telnet  program,  which  allowed  us  our  first 
access  to  the  ARPANET,  and  the  development  of  other  supporting  software. 

In  July,  with  cooperation  from  the  staff  at  the  1 1.  L I AC-  IV  Center  for 
Advanced  Computation,  our  first  running  ELF  TELNET  server  was  installed.  It  is 
a copy  of  the  ELF  software,  custom  tailored  to  our  system.  The  current  version 
supports  TELNET  through  our  terminals  and  FTP  to  our  line-printer.  Since  we 
expected  to  replace  the  system  with  NETWORK-UN I X,  we  did  not  expand  the  ELF  to 
use  the  other  devices  in  our  system.  Since  UNIX  did  not  have  the  network  soft- 
ware installed  when  it  first  arrived,  ELF  served  as  our  primary  entry  to  the 
network.  However,  when  the  full,  NETWORK-UNIX  system  was  operational,  ELF 
became  obsolete,  and  is  rearly  used,  except  in  those  instances  where  hardware 
problems  disallow  the  use  of  UNIX,  yet  still  allow  ELF  to  run. 

Our  use  of  the  ARPANET  remains  small.  Most  of  our  research  has  been 
carried  out  either  on  our  own  PDP  system,  or  at  the  Laboratory  for  the  Applica- 
tion of  Remote  Sensing  and  the  Purdue  campus  central  facility.  These  latter 
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sites  have  the  large,  scientific  machines  required  for  our  research  (two  CDC 

SSO^'s  and  an  IBM  360).  In  the  future,  we  plan  to  use  'he  360/91  at  UCLA  and 
the  I LL I AC- IV. 

The  development  of  PDP-11  software  over  the  year  can  be  considered 

| 

quarter  by  qua* ter.  During  the  first  quarter  our  primary  efforts  were  directed 
towards  checking  out  the  hardware  and  getting  the  DOS  system  operational.  This 
was  complicated  by  the  fact  that  there  were  many  hardware  problems,  and  by  the 
fact  that  the  DOS  system  did  not  arrive  until  the  middle  of  January.  This  de- 
layed the  correction  of  the  hardware  problems  because  many  of  them  were  not 
discovered  until  after  DOS  was  delivered.  We  also  spent  time  the  first  quarter 
familiarizing  ourselves  with  the  operation  of  the  equipment  and  the  DOS  soft- 
ware. During  the  serond  quarter  we  concerned  ourselves  with  developing  basic 
software  for  DOS  which  allowed  access  to  the  non-standard  devices,  like  the 
COMTAL  color  display  and  the  VERSATEK  electrostatic  printer,  and  other  basic 
software  to  use  these  devices.  We  also  were  investigating  the  various  network 
software  products  which  were  available  at  that  time.  We  were  considering  both 
ANTS  and  ELF.  By  the  end  of  the  second  quarter,  we  had  enough  information  about 
ANTS  and  ELF,  and  enough  experience  with  the  ARPANET  to  be  able  to  begin  working 
on  the  installation  of  a network  software  system.  Most  of  the  third  quarter 
was  spent  first  by  debugging  the  IMP-interface  hardware,  and  working  on  getting 
ELF  running  on  our  machine.  The  software  portion  of  this  work  was  complicated 
by  the  fact  that  we  did  rj1:  have  access  to  the  network  except  via  a noisy  line 
to  Illinois.  Our  original  plan  was  to  get  a magnetic  tape  containing  the  E'.F 
sources,  and  compile  them  on  our  PDP-11,  using  DOS.  Unfortunately,  the  DOS 
system  was  not  able  to  compile  the  (very  large)  ELF  sources,  and  this  was  not 
discovered  until  several  abortive  attempts  were  made.  About  the  time  that  the 
hardware  was  finished  (late  May),  we  accepted  an  offer  from  the  ILLIAC-IV  IAC 
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to  help  us  (jet  a copy  of  the  f r ELF  running  on  our  machine.  With  their  assis- 
tance, we  got  the  network  software  running  in  early  July.  It  was  about  this 
time  when  we  heard  about  UNIX,  and  the  network  software  developed  for  it  by 
University  of  Illinois.  We  researched  it,  and  decided  that  for  a machine  as 
large  as  ours,  that  ELF  wasted  most  of  the  equipment  we  had.  UNIX,  on  the 
other  hand,  allowed  us  to  access  the  network,  and  to  perform  local  processing, 
simultaneously.  Thus  most  of  the  last  quarter  was  spent  installing,  improving 
and  testing  the  UNIX  system.  We  were  very  pleased.  UNIX,  used  only  for  net- 
work access  turned  out  to  be  better  than  ELF,  and  used  only  for  local  process- 
ing was  better  than  DOS  or  RSX.  Taken  as  a whole,  it  seems  to  be  the  best 
system  available  for  a moderately  sized  system  like  ours.  This  opinion  seem;, 
to  oe  shared  by  otht  PDP-11  network  users,  as  several  of  the  old  ELF-fol  lov/ers 
have  switched  to  NETWORK-UNI  X.  UNIX  supports  FORTRAN,  PDP-11  assemply  language, 
and  a frgh-level  language,  similar  to  PL/1,  called  "C".  It  also  supports  several 
"compiler-compilers"  and  has  a computer  language  translator,  which  will  be  a 
tremendous  aid  to  us  as  we  translate  out  CDC  and  IBM  software  to  run  on  the 
PDP-11.  That  UNIX  is  an  easy  system  to  use  is  supported  by  the  fact  that  to 
develop  the  software  to  run  our  COMTAL  display  and  the  electrostatic  printer 
has  only  about  20%  of  the  development  time  for  the  equivalent  software  for  DOS. 

This  next  year  we  are  continuing  the  development  of  more  UNIX  support 
software,  and  converting  old  programs  to  run  under  UNIX;  and  in  particular,  we 
are  starting  to  do  most  of  our  new  research  projects  on  the  PDP-11  UNIX  system. 
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Mnnu  f oc  ttircr 


Beehive  Elect. 
Tex.  Inst. 
Digi-Data 

DEC 

DEC 

Fabr i tek 

Versatek 

Comtal 

Data  Printer 
T rue-Data 
Tektron i x 
DEC 


Current  Equipment  Configuration 
February,  1976 

Descr Ipt Ion 
"Super-Bee"  Terminals 

"Silent  700"  Terminals 

Industry  standard  magnetic  tape  system; 

2,  9"track  and  1,  7-track  drives;  one  each 
NRZI  and  phase-encoded  formatters/controllers 

Dual-drive  DECtape  unit 

RP03  disk  drive  (40  million  characters) 

96K-word  auxiliary  memory  system 
(6^K  bought  by  ARPA,  32K  by  NASA) 

Electrostatic  matrix  printer 

Color  picture  display 

132  column,  600  L.P.M.  line  printer 

Punched  card  reader 

Model  4010,  graphics  display 

PDP-11/45  computer  system;  system  Includes: 
32K  memory 

FPP-11  floating  point  processor  (NSF  money) 
H960  extension  mounting  cabinet 
3 - small  peripheral  mountings  blocks  (DD-11) 

1 UN! BUS  repeater/expander 

DH11,  16-line  terminal  multiplexor 

KWll-p  programmed  clock 

"ANTS"  - type  PDP-11/ IMP  interface 
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